### PREDICTING ELECTION OCCURENCE FROM ECONOMIC CRISES -- ANNUAL DATA
require(blme)
cat("Figure 9: Endogenous election timing: will take a while.")
load("Data/annualdata.Rdata")
load("Data/nelda.Rdata")

nelda.pres <- nelda[which(nelda$legelec==0 & nelda$round==1), c("country", "year")]
nelda.pres$presidential <- 1
qog <- merge(qog, nelda.pres, by=c("country", "year"), all.x=TRUE)
qog$presidential <- 1*(!is.na(qog$presidential))

nelda.parl <- nelda[which(nelda$legelec==1 & nelda$round==1), c("country", "year")]
nelda.parl$parliamentary <- 1
qog <- merge(qog, nelda.parl, by=c("country", "year"), all.x=TRUE)
qog$parliamentary <- 1*(!is.na(qog$parliamentary))
qog <- qog[order(qog$country, qog$year), ]
lcumsum <- function(x) cumsum(c(0, x[-length(x)]))
qog$presspell <- unlist(with(qog, tapply(presidential, country, lcumsum)))
qog$parlspell <- unlist(with(qog, tapply(parliamentary, country, lcumsum)))

qog$parl.time <- unlist(with(qog, tapply(parliamentary, paste(country, parlspell), function(x) 1:length(x))))
qog$pres.time <- unlist(with(qog, tapply(presidential, paste(country, presspell), function(x) 1:length(x))))

qog$presels <- with(qog, tapply(presidential,  paste(country, presspell), sum)[paste(country, presspell)])

qog$parlels <- with(qog, tapply(parliamentary,  paste(country, parlspell), sum)[paste(country, parlspell)])

qog <- subset(qog, agedem_t1 <= 5 & !is.na(growthc_t1))

L <- 3
sfit1 <- bglmer(presidential ~ ns(growthc_t1, L) + (1 + ns(pres.time, L)|country), data=qog, family=binomial(link="probit"), glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+07)), cov.prior = invwishart(df = 5, scale = diag(0.5, L+1)))
sfit2 <- update(sfit1, subset=presels > 0)         
sfit3 <- update(sfit1, parliamentary ~ ns(growthc_t1, L) + (1 + ns(parl.time, L)|country))
sfit4 <- update(sfit3, subset=parlels > 0)

pplot <- function(fit, title=NULL){

  q <- qog$growthc_t1
  x <- seq(min(q), max(q), length=100)
  b <- fixef(fit)
  X <- cbind(1, predict(ns(qog$growthc_t1, L), x))

  bhat <- rmnorm(1e+4, b, as.matrix(vcov(fit)))
  p <- pnorm(t(X%*%t(bhat)))
  P <- t(apply(p,2,function(x) quantile(x, c(0.05,0.5, 0.95))))

  data <- data.frame(x=x, P)
  data <- melt(data, id="x")

  p <- ggplot(data=data, aes(x=x, y=value, linetype=variable)) + geom_line(size=1.5, colour="black") + scale_linetype_manual(values=c(3,1,3)) + labs(x = "Lagged pc GDP growth (standardized)", y="Predicted probability of elections")  + scale_colour_manual(values=c("black", "black")) +  theme_bw() + theme(
  legend.position = "none", 
  panel.grid.minor = element_blank(),
  panel.background = element_rect(colour = "white"),
  axis.text.x = element_text(size=12),
  axis.text.y = element_text(size=12),
  plot.title=element_text(face="bold"),
  panel.grid.minor = element_blank(),
  panel.background = element_rect(colour = "white")
) + labs(title=title)
  p + scale_y_continuous(breaks = seq(0,.5,by=.1), limits = c(0, .55))
}

pdf("../graphs/growthelections_paper.pdf", width=10, height=5)
grid.arrange(pplot(sfit1, "Presidential elections"), pplot(sfit3, "Parliamentary elections"), ncol=2)
dev.off()

pdf("../graphs/growthelections.pdf", width=8, height=8)
grid.arrange(pplot(sfit1, "Presidential (all sample)"), pplot(sfit2, "Parliamentary (all sample)"), pplot(sfit3, "Presidential (excluding countries that \never held presidential elections)"), pplot(sfit4, "Parliamentary (excluding countries that \never held parliamentary elections)"), ncol=2)
dev.off()